##### ##################################################################
#####
#####   Input: Clean experimental data
#####   Output: Analysis and figures for the supplementary materials
#####

setwd("/Users/lotte/Dropbox/projects/productivity_experiment/replication")

rm(list = ls())

library(data.table) # CRAN v1.14.2
library(plyr) # CRAN v1.8.6
library(dplyr) # CRAN v1.0.9
library(tidyverse) # CRAN v1.3.1
library(cjoint) # CRAN v2.1.0
library(ggplot2) # CRAN v3.3.6
library(broom) # CRAN v0.7.11
library(patchwork) # CRAN v1.1.1
library(estimatr) # CRAN v0.30.6

# Load data

load("data/productivity.Rdata")

# Unconditional - all attributes

preference_bivariate <-
  lm(
    preference ~ mp_gender + mp_party + committee_membership +
      issue_campaigning + voting_legislation + constituency_responsiveness,
    data = productivity
  )

performance_bivariate <-
  lm(
    perceived_performance ~ mp_gender + mp_party + committee_membership +
      issue_campaigning + voting_legislation + constituency_responsiveness,
    data = productivity
  )

electability_bivariate <-
  lm(
    electability ~ mp_gender + mp_party + committee_membership +
      issue_campaigning + voting_legislation + constituency_responsiveness,
    data = productivity
  )

save(preference_bivariate,
     performance_bivariate,
     electability_bivariate,
     file = "tables/unconditional.Rdata")

# Conditional - MP gender

## Continuous objective productivity

electability <-
  lm(electability ~ mp_gender * objective_performance, data = productivity)
job_performance <-
  lm(perceived_performance ~ mp_gender * objective_performance,
     data = productivity)
preference <-
  lm(preference ~ mp_gender * objective_performance, data = productivity)

## Binary productivity attributes

electability_levels <-
  lm(
    electability ~ mp_gender * committee_membership + mp_gender * issue_campaigning + mp_gender *
      voting_legislation +
      mp_gender * constituency_responsiveness,
    data = productivity
  )

job_performance_levels <-
  lm(
    perceived_performance ~ mp_gender * committee_membership + mp_gender * issue_campaigning + mp_gender *
      voting_legislation +
      mp_gender * constituency_responsiveness,
    data = productivity
  )

preference_levels <-
  lm(
    preference ~ mp_gender * committee_membership + mp_gender * issue_campaigning + mp_gender *
      voting_legislation +
      mp_gender * constituency_responsiveness,
    data = productivity
  )

save(
  electability,
  job_performance,
  preference,
  electability_levels,
  job_performance_levels,
  preference_levels,
  file = "tables/conditional_mp_gender.Rdata"
)

## Non-linear

job_performance_0 <-
  lm(perceived_performance ~ mp_gender * performance_0_dummy, data = productivity)
job_performance_1 <-
  lm(perceived_performance ~ mp_gender * performance_1_dummy, data = productivity)
job_performance_2 <-
  lm(perceived_performance ~ mp_gender * performance_2_dummy, data = productivity)
job_performance_3 <-
  lm(perceived_performance ~ mp_gender * performance_3_dummy, data = productivity)
job_performance_4 <-
  lm(perceived_performance ~ mp_gender * performance_4_dummy, data = productivity)
electability_0 <-
  lm(electability ~ mp_gender * performance_0_dummy, data = productivity)
electability_1 <-
  lm(electability ~ mp_gender * performance_1_dummy, data = productivity)
electability_2 <-
  lm(electability ~ mp_gender * performance_2_dummy, data = productivity)
electability_3 <-
  lm(electability ~ mp_gender * performance_3_dummy, data = productivity)
electability_4 <-
  lm(electability ~ mp_gender * performance_4_dummy, data = productivity)
preference_0 <-
  lm(preference ~ mp_gender * performance_0_dummy, data = productivity)
preference_1 <-
  lm(preference ~ mp_gender * performance_1_dummy, data = productivity)
preference_2 <-
  lm(preference ~ mp_gender * performance_2_dummy, data = productivity)
preference_3 <-
  lm(preference ~ mp_gender * performance_3_dummy, data = productivity)
preference_4 <-
  lm(preference ~ mp_gender * performance_4_dummy, data = productivity)


save(
  job_performance_0,
  job_performance_1,
  job_performance_2,
  job_performance_3,
  job_performance_4,
  electability_0,
  electability_1,
  electability_2,
  electability_3,
  electability_4,
  preference_0,
  preference_1,
  preference_2,
  preference_3,
  preference_4,
  file = "tables/conditional_nonlinear.Rdata"
)

## Conditional and unconditional - party

electability_unconditional <-
  lm(electability ~ party_congruence, data = productivity)
electability <-
  lm(electability ~ party_congruence * objective_performance, data = productivity)
job_performance_unconditional <-
  lm(perceived_performance ~ party_congruence, data = productivity)
job_performance <-
  lm(perceived_performance ~ party_congruence * objective_performance,
     data = productivity)
preference_unconditional <-
  lm(preference ~ party_congruence, data = productivity)
preference <-
  lm(preference ~ party_congruence * objective_performance, data = productivity)

electability_levels <-
  lm(
    electability ~ party_congruence * committee_membership + party_congruence *
      issue_campaigning +
      party_congruence * voting_legislation + party_congruence *
      constituency_responsiveness,
    data = productivity
  )

job_performance_levels <-
  lm(
    perceived_performance ~ party_congruence * committee_membership + party_congruence *
      issue_campaigning +
      party_congruence * voting_legislation + party_congruence *
      constituency_responsiveness,
    data = productivity
  )

preference_levels <-
  lm(
    preference ~ party_congruence * committee_membership + party_congruence *
      issue_campaigning +
      party_congruence * voting_legislation + party_congruence *
      constituency_responsiveness,
    data = productivity
  )

save(
  electability,
  job_performance,
  preference,
  electability_levels,
  job_performance_levels,
  preference_levels,
  electability_unconditional,
  job_performance_unconditional,
  preference_unconditional,
  file = "tables/conditional_party_congruence.Rdata"
)

# Robustness

## Mixed-sex pairings only - AMCE

male_female <-
  productivity[productivity$gender_pairing == "Man-Woman", ]

results_male_female <-
  amce(
    preference ~ mp_gender + mp_party + committee_membership +
      issue_campaigning + voting_legislation + constituency_responsiveness,
    respondent.id = NULL,
    data = male_female
  )

pdf("plots/figure_S1_amce_male_female.pdf", 7, 4)
plot(
  results_male_female,
  xlab = "Change in Pr(Preferred Member of Parliament)",
  attribute.names = c(
    "Committee Membership",
    "Constituency Responsiveness",
    "Issue Campaigning",
    "MP Gender",
    "MP Party",
    "Voting and Legislation"
  ),
  main = "",
  text.size = 10
)
dev.off()

## Mixed-sex pairings only - Continuous

model_1 <-
  lm_robust(
    preference ~ objective_performance * mp_gender,
    data = male_female,
    clusters = ID,
    se_type = "stata"
  )
model_2 <-
  lm_robust(
    electability ~ objective_performance * mp_gender,
    data = male_female,
    clusters = ID,
    se_type = "stata"
  )
model_3 <-
  lm_robust(
    perceived_performance ~ objective_performance * mp_gender,
    data = male_female,
    clusters = ID,
    se_type = "stata"
  )

new_data_man <- data.frame(mp_gender = "Man",
                           objective_performance = c(0:4))
new_data_woman <- data.frame(mp_gender = "Woman",
                             objective_performance = c(0:4))

model_1_man <-
  predict(model_1, interval = "confidence", newdata = new_data_man)
model_1_woman <-
  predict(model_1, interval = "confidence", newdata = new_data_woman)
model_2_man <-
  predict(model_2, interval = "confidence", newdata = new_data_man)
model_2_woman <-
  predict(model_2, interval = "confidence", newdata = new_data_woman)
model_3_man <-
  predict(model_3, interval = "confidence", newdata = new_data_man)
model_3_woman <-
  predict(model_3, interval = "confidence", newdata = new_data_woman)

model_1_for_plot <-
  as.data.frame(rbind(model_1_man$fit, model_1_woman$fit)) %>%
  mutate(
    gender = c(
      "Man",
      "Man",
      "Man",
      "Man",
      "Man",
      "Woman",
      "Woman",
      "Woman",
      "Woman",
      "Woman"
    ),
    outcome = c("MP preference"),
    objective_performance = c(0, 1, 2, 3, 4, 0, 1, 2, 3, 4)
  )

model_2_for_plot <-
  as.data.frame(rbind(model_2_man$fit, model_2_woman$fit)) %>%
  mutate(
    gender = c(
      "Man",
      "Man",
      "Man",
      "Man",
      "Man",
      "Woman",
      "Woman",
      "Woman",
      "Woman",
      "Woman"
    ),
    outcome = c("Electability"),
    objective_performance = c(0, 1, 2, 3, 4, 0, 1, 2, 3, 4)
  )

model_3_for_plot <-
  as.data.frame(rbind(model_3_man$fit, model_3_woman$fit)) %>%
  mutate(
    gender = c(
      "Man",
      "Man",
      "Man",
      "Man",
      "Man",
      "Woman",
      "Woman",
      "Woman",
      "Woman",
      "Woman"
    ),
    outcome = c("Job Performance"),
    objective_performance = c(0, 1, 2, 3, 4, 0, 1, 2, 3, 4)
  )

models_for_plot <- rbind(model_1_for_plot,
                         model_2_for_plot,
                         model_3_for_plot)

male_female_only <-
  ggplot(models_for_plot,
         aes(y = fit, x = objective_performance, group = gender)) +
  geom_ribbon(aes(
    ymin = lwr,
    ymax = upr,
    linetype = gender,
    fill = gender
  ), alpha = 0.3) +
  theme_bw() +
  scale_fill_manual(values = c("lightseagreen", "grey40")) +
  scale_colour_manual(values = c("lightseagreen", "grey40")) +
  geom_line(aes(
    group = gender,
    color = gender,
    linetype = gender
  )) +
  labs(color = "MP gender",
       linetype = "MP gender",
       fill = "MP gender") +
  facet_wrap( ~ outcome, scales = "free") +
  ylab("") + xlab("Objective performance") +
  theme(
    axis.text = element_text(size = 9),
    strip.text.x = element_text(size = 9),
    legend.title = element_text(size = 10),
    axis.title = element_text(size = 9)
  )
ggsave(
  "plots/figure_S2_interactions_male_female_pairings.pdf",
  male_female_only,
  width = 10,
  height = 5
)

# HTEs

## Gender congruence
### Marginal means

library(cregg) # CRAN v0.4.0

mm_gender_congruence <-
  cj(
    productivity,
    preference ~ mp_gender + mp_party + committee_membership +
      issue_campaigning + voting_legislation + constituency_responsiveness,
    respondent.id = "ID",
    estimate = "mm",
    by =  ~ gender_congruence
  )

pdf("plots/figure_S3_mm_gender_congruence.pdf", 7, 5)
plot(
  mm_gender_congruence,
  group = "gender_congruence",
  vline = 0.5,
  text.size = 9
)
dev.off()

### Continuous

model_6 <-
  lm_robust(
    preference ~ objective_performance * gender_congruence,
    data = productivity,
    clusters = ID,
    se_type = "stata"
  )
model_7 <-
  lm_robust(
    perceived_performance ~ objective_performance * gender_congruence,
    data = productivity,
    clusters = ID,
    se_type = "stata"
  )
model_8 <-
  lm_robust(
    electability ~ objective_performance * gender_congruence,
    data = productivity,
    clusters = ID,
    se_type = "stata"
  )

new_data_congruent <- data.frame(gender_congruence = "TRUE",
                                 objective_performance = c(0:4))
new_data_incongruent <- data.frame(gender_congruence = "FALSE",
                                   objective_performance = c(0:4))
model_6_congruent <-
  predict(model_6, interval = "confidence", newdata = new_data_congruent)
model_6_incongruent <-
  predict(model_6, interval = "confidence", newdata = new_data_incongruent)
model_6_congruent_for_plot <-
  as.data.frame(rbind(model_6_congruent$fit, model_6_incongruent$fit)) %>%
  mutate(
    congruence = c(
      "Congruent",
      "Congruent",
      "Congruent",
      "Congruent",
      "Congruent",
      "Incongruent",
      "Incongruent",
      "Incongruent",
      "Incongruent",
      "Incongruent"
    ),
    outcome = c("MP Preference"),
    objective_performance = c(0, 1, 2, 3, 4, 0, 1, 2, 3, 4)
  )
model_7_congruent <-
  predict(model_7, interval = "confidence", newdata = new_data_congruent)
model_7_incongruent <-
  predict(model_7, interval = "confidence", newdata = new_data_incongruent)
model_7_congruent_for_plot <-
  as.data.frame(rbind(model_7_congruent$fit, model_7_incongruent$fit)) %>%
  mutate(
    congruence = c(
      "Congruent",
      "Congruent",
      "Congruent",
      "Congruent",
      "Congruent",
      "Incongruent",
      "Incongruent",
      "Incongruent",
      "Incongruent",
      "Incongruent"
    ),
    outcome = c("Job Performance"),
    objective_performance = c(0, 1, 2, 3, 4, 0, 1, 2, 3, 4)
  )

model_8_congruent <-
  predict(model_8, interval = "confidence", newdata = new_data_congruent)
model_8_incongruent <-
  predict(model_8, interval = "confidence", newdata = new_data_incongruent)
model_8_congruent_for_plot <-
  as.data.frame(rbind(model_8_congruent$fit, model_8_incongruent$fit)) %>%
  mutate(
    congruence = c(
      "Congruent",
      "Congruent",
      "Congruent",
      "Congruent",
      "Congruent",
      "Incongruent",
      "Incongruent",
      "Incongruent",
      "Incongruent",
      "Incongruent"
    ),
    outcome = c("Electability"),
    objective_performance = c(0, 1, 2, 3, 4, 0, 1, 2, 3, 4)
  )

gender_congruence <-
  rbind(
    model_6_congruent_for_plot ,
    model_7_congruent_for_plot,
    model_8_congruent_for_plot
  )

gender_congruence_plot <-
  ggplot(gender_congruence,
         aes(y = fit, x = objective_performance, group = congruence)) +
  geom_ribbon(aes(
    ymin = lwr,
    ymax = upr,
    linetype = congruence,
    fill = congruence
  ),
  alpha = 0.3) +
  theme_bw() +
  scale_fill_manual(values = c("lightseagreen", "grey40")) +
  scale_colour_manual(values = c("lightseagreen", "grey40")) +
  geom_line(aes(
    group = congruence,
    color = congruence,
    linetype = congruence
  )) +
  labs(color = "Gender congruence",
       linetype = "Gender congruence",
       fill = "Gender congruence") +
  facet_wrap( ~ outcome, scales = "free") +
  ylab("") + xlab("Objective performance") +
  theme(
    axis.text = element_text(size = 9),
    strip.text.x = element_text(size = 9),
    legend.title = element_text(size = 9),
    axis.title = element_text(size = 9)
  )
ggsave(
  "plots/figure_S5_gender_congruence.pdf",
  gender_congruence_plot,
  width = 10,
  height = 5
)

## Respondent gender
### Marginal means

mm_gender <-
  cj(
    productivity,
    preference ~ mp_gender + mp_party + committee_membership +
      issue_campaigning + voting_legislation + constituency_responsiveness,
    respondent.id = "ID",
    estimate = "mm",
    by =  ~ respondent_gender
  )

pdf("plots/figure_S4_mm_gender.pdf", 7, 5)
plot(mm_gender,
     group = "respondent_gender",
     vline = 0.5,
     text.size = 9)
dev.off()

### Continuous

men_only <- productivity[productivity$respondent_gender == "Man", ]
women_only <- productivity[productivity$respondent_gender == "Woman", ]

model_1_men <-
  lm_robust(
    perceived_performance ~ objective_performance * mp_gender,
    data = men_only,
    clusters = ID,
    se_type = "stata"
  )
model_2_men <-
  lm_robust(
    preference ~ objective_performance * mp_gender,
    data = men_only,
    clusters = ID,
    se_type = "stata"
  )
model_3_men <-
  lm_robust(
    electability ~ objective_performance * mp_gender,
    data = men_only,
    clusters = ID,
    se_type = "stata"
  )

new_data_man <- data.frame(mp_gender = "Man",
                           objective_performance = c(0:4))
new_data_woman <- data.frame(mp_gender = "Woman",
                             objective_performance = c(0:4))

model_1_man <-
  predict(model_1_men, interval = "confidence", newdata = new_data_man)
model_1_woman <-
  predict(model_1_men, interval = "confidence", newdata = new_data_woman)
model_2_man <-
  predict(model_2_men, interval = "confidence", newdata = new_data_man)
model_2_woman <-
  predict(model_2_men, interval = "confidence", newdata = new_data_woman)
model_3_man <-
  predict(model_3_men, interval = "confidence", newdata = new_data_man)
model_3_woman <-
  predict(model_3_men, interval = "confidence", newdata = new_data_woman)

model_1_for_plot <-
  as.data.frame(rbind(model_1_man$fit, model_1_woman$fit)) %>%
  mutate(
    gender = c(
      "Man",
      "Man",
      "Man",
      "Man",
      "Man",
      "Woman",
      "Woman",
      "Woman",
      "Woman",
      "Woman"
    ),
    outcome = c("Job Performance"),
    voter_gender = c("Man Voter"),
    objective_performance = c(0, 1, 2, 3, 4, 0, 1, 2, 3, 4)
  )

model_2_for_plot <-
  as.data.frame(rbind(model_2_man$fit, model_2_woman$fit)) %>%
  mutate(
    gender = c(
      "Man",
      "Man",
      "Man",
      "Man",
      "Man",
      "Woman",
      "Woman",
      "Woman",
      "Woman",
      "Woman"
    ),
    outcome = c("MP Preference"),
    voter_gender = c("Man Voter"),
    objective_performance = c(0, 1, 2, 3, 4, 0, 1, 2, 3, 4)
  )

model_3_for_plot <-
  as.data.frame(rbind(model_3_man$fit, model_3_woman$fit)) %>%
  mutate(
    gender = c(
      "Man",
      "Man",
      "Man",
      "Man",
      "Man",
      "Woman",
      "Woman",
      "Woman",
      "Woman",
      "Woman"
    ),
    outcome = c("Electability"),
    voter_gender = c("Man Voter"),
    objective_performance = c(0, 1, 2, 3, 4, 0, 1, 2, 3, 4)
  )

model_1_women <-
  lm_robust(
    perceived_performance ~ objective_performance * mp_gender,
    data = women_only,
    clusters = ID,
    se_type = "stata"
  )
model_2_women <-
  lm_robust(
    preference ~ objective_performance * mp_gender,
    data = women_only,
    clusters = ID,
    se_type = "stata"
  )
model_3_women <-
  lm_robust(
    electability ~ objective_performance * mp_gender,
    data = women_only,
    clusters = ID,
    se_type = "stata"
  )

model_1_man_w <-
  predict(model_1_women, interval = "confidence", newdata = new_data_man)
model_1_woman_w <-
  predict(model_1_women, interval = "confidence", newdata = new_data_woman)
model_2_man_w <-
  predict(model_2_women, interval = "confidence", newdata = new_data_man)
model_2_woman_w <-
  predict(model_2_women, interval = "confidence", newdata = new_data_woman)
model_3_man_w <-
  predict(model_3_women, interval = "confidence", newdata = new_data_man)
model_3_woman_w <-
  predict(model_3_women, interval = "confidence", newdata = new_data_woman)

model_1_for_plot_women <-
  as.data.frame(rbind(model_1_man_w$fit, model_1_woman_w$fit)) %>%
  mutate(
    gender = c(
      "Man",
      "Man",
      "Man",
      "Man",
      "Man",
      "Woman",
      "Woman",
      "Woman",
      "Woman",
      "Woman"
    ),
    outcome = c("Job Performance"),
    voter_gender = c("Woman Voter"),
    objective_performance = c(0, 1, 2, 3, 4, 0, 1, 2, 3, 4)
  )

model_2_for_plot_women  <-
  as.data.frame(rbind(model_2_man_w$fit, model_2_woman_w$fit)) %>%
  mutate(
    gender = c(
      "Man",
      "Man",
      "Man",
      "Man",
      "Man",
      "Woman",
      "Woman",
      "Woman",
      "Woman",
      "Woman"
    ),
    outcome = c("MP Preference"),
    voter_gender = c("Woman Voter"),
    objective_performance = c(0, 1, 2, 3, 4, 0, 1, 2, 3, 4)
  )

model_3_for_plot_women  <-
  as.data.frame(rbind(model_3_man_w$fit, model_3_woman_w$fit)) %>%
  mutate(
    gender = c(
      "Man",
      "Man",
      "Man",
      "Man",
      "Man",
      "Woman",
      "Woman",
      "Woman",
      "Woman",
      "Woman"
    ),
    outcome = c("Electability"),
    voter_gender = c("Woman Voter"),
    objective_performance = c(0, 1, 2, 3, 4, 0, 1, 2, 3, 4)
  )

respondent_gender_for_plot <- rbind(
  model_1_for_plot,
  model_2_for_plot,
  model_3_for_plot,
  model_1_for_plot_women,
  model_2_for_plot_women,
  model_3_for_plot_women
)

hte_gender <-
  ggplot(respondent_gender_for_plot,
         aes(y = fit, x = objective_performance, group = gender)) +
  geom_ribbon(aes(
    ymin = lwr,
    ymax = upr,
    linetype = gender,
    fill = gender
  ), alpha = 0.3) +
  theme_bw() +
  scale_fill_manual(values = c("lightseagreen", "grey40")) +
  scale_colour_manual(values = c("lightseagreen", "grey40")) +
  geom_line(aes(
    group = gender,
    color = gender,
    linetype = gender
  )) +
  labs(color = "MP gender",
       linetype = "MP gender",
       fill = "MP gender") +
  facet_grid(voter_gender ~ outcome, scales = "free") +
  ylab("") + xlab("Objective performance") +
  theme(
    axis.text = element_text(size = 9),
    strip.text.x = element_text(size = 9),
    legend.title = element_text(size = 10),
    axis.title = element_text(size = 9)
  )
ggsave("plots/figure_S6_hte_gender.pdf",
       hte_gender,
       width = 10,
       height = 5)

